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ABSTRACT 

We study the migration of solid bodies in turbulent protoplanetary accretion discs by 
means of global MHD simulations. The bodies range in size from 5 centimetres up 
to 1 metre, and so include objects whose migration is expected to be the most rapid 
due to gas drag interaction with the disc. As they drift inward through the disc, some 
of them are trapped in regions where gas pressure maxima are created by long lived 
anticyclonic vortices. This accumulation is very efficient, locally increasing the dust- 
to-gas ratio by a factor > 100 in some cases. We discuss the possible implications of 
this result for theories of planet formation. 
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1 INTRODUCTION 

A crucial element of any planet formation theory is under- 
standing how the dust component of protoplanetary discs 
builds up to form larger objects that eventually become 
planetesimals and planets. This has obvious relevance to the 
formation of rocky terrestrial planets and asteroids, and ac- 
cording to the core-accretion model of gas giant formation 
(e.g. Pollack et al. 1996), is a key stage in the formation of 
giant planets. 

The growth and evolution of solid bodies of all sizes 
is therefore an important topic to study. A key issue re- 
lates to the growth and survivability of metre-sized objects, 
which are ex pected to undergo rapi d inward migration due 
to gas drag iWeidenschilling||l977l) . Such bodies are pre- 
dicted to migrate into the central star within ~ 100 years 
of formation, raising questions about how planetary systems 
form at all. Models of planet formation in smooth, laminar 
protoplanetary discs, beginning with small dust grains that 
must grow through this size range, require a radially ex- 
tend ed reservoir of solid matter to overcome this problem 
(e.g. IWeidenschilling fc Cuzzilll993l) . 

In this letter, we investigate the effect MHD turbulence 
resulting from the magneto-rotational instability (MRI) in 
the disc might have on the dynamics of such objects. It seems 
likely that the MRI is the source of turbulence which trans- 
ports angular momentum in protoplanetary discs. The ve- 
locity and density fluctuations that result will modify the 
drag force exerted on solid bodies and will affect their dy- 
namics. These fluctuations can range from being relatively 
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small variations, up to being substantial modifications of 
the disc radial profil e caused by radial variations in the tur- 
bulen t stresses (e.g. lHawlevll200ft ISteinacker fc Papaloizoul 
12002ft . 

Because of the global nature of this problem, in which 
solid bodies will drift through a substantial fraction of the 
disc r adius, we set up global proto planetary disc simulations. 
As in Papaloizou fc Nelson! (: 2003). we have developed cylin- 
drical discs models to limit the large computational cost 
of these calculations. We not e that there has been a re- 
cent s tudy of this problem bv lJohansen. Klahr fc Hennind 
(2005), who find transient trapping of solid bodies in den- 
sity maxima generated by MHD turbulence. The approach 
used in their study, however, was local rather than the global 
approach we have taken. 

The plan of the paper is as follows. In section 2, we 
describe our numerical schemes, and in section 3 we describe 
our results. We discuss their consequences and limitations in 
section 4. 



2 SIMULATIONS 
2.1 Algorithms 

Two similar MHD Eulerian codes were used to com- 
pute the simulat i ons p resented in this p aper: GLO BAL 
ilHawlev fc Stond 119951) and NIRVANA iZierier fc Yorkd 
Il997ft . Both use finite difference techniques combined with 
the Constrained Transport algorithm to evolve the MHD 
equations. They differ, however, in their treatment of the 
solid phase: GLOBAL treats the solid bodies as a second 
fluid; NIRVANA treats the solid bodies as particles. This 
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reflects the sizes of the objects considered in each of our 
simulations. 



2.1.1 The two fluid approach 

GLOBAL was extended to describe the dust particle evolu- 
tion as a second, pressureless fluid. This second fluid inter- 
acts with the gas via the drag force exerted by the latter. 



The drag force Fd rag takes the form 



Fdr 



— (Vg - V d ) 



(1) 



Here m p is the mass of the particle, (v g — Vd) is the gas ve- 
locity relative to the solid component, and r s is the dust 
stopping time. The small dust particles examined using 
GLOBAL are in the Epstein regime, so t s is defined by 

pC s 

where p is the gas density, c s the speed of sound, p s the 
density of the dust particles and a their size. Because of 
the short stopping time of the small particles studied in this 
work, the effect of the drag force is computed implicitly. 



Model 


dust 


Number of 


Resolution 


Particle size 




description 


dust particles 




(in cm) 


Gl 


fluid 


No particle 


(260, 152, 44) 


5 


G2 


fluid 


No particle 


(260, 152, 44) 


25 


Nl 


N-body 


3000 


(260, 604, 44) 


100 


N2 


N-body 


500 


(260, 604, 44) 


100 



Table 1. Model parameters of the runs presented here: Column 1 
gives the name of the model. The first two (Gl and G2) are per- 
formed with GLOBAL and the last two (Nl and N2) with NIR- 
VANA. The algorithm used to evolve the dust component (fluid 
vs N— body) is described in column 2. For the particle approach, 
column 3 gives the number of dust particles used. Finally, col- 
umn 4 gives the resolution of the run and column 5 the size of 
the solid bodies, in centimetres. 




2.1.2 The N-body approach 

The two fluid approximation breaks down for large solid 
bodies, so we adapted NIRVANA to include coupling be- 
tween the fluid and individual solid objects represented by 
particles. We consider metre-sized bodies, whose drag force 
interaction is given by Stokes Law for the disc parameters 
we have adopted. The stopping time, r s , then becomes (e.g. 
IWeidenschiliinslll977l) 



_ 2p s a 

where r\ is the gas viscosity defined by 



(3) 



(4) 



and the mean free path of molecules A = (nH 2 o"H 2 ) _1 - We 
adopt a value of <th 2 = 10 -15 cm 2 for the collision cross sec- 
tion of molecular hydrogen, and assume that nu 2 = p/m^, 
wher e is the mass of a hydrogen molecule (e.g. [Rafikov 
12004) . We employ linear interpolation, using information 
from the surrounding eight cells, to obtain relevant phys- 
ical quantities at the location of the particles. 

In one simulation performed (model Nl), the parti- 
cles interacted with the disc through the drag force only. 
A second simulation (model N2) was performed in which 
the particles also experienced gravitational interaction with 
the disc. In this case the axisymmetric component of disc 
self-gravity was included to ensure that the velocity differ- 
ence betw een solid particles and the gas was c omputed cor- 
rectly (see iFromang. Terauem fc Nelsonll2005l for a descrip- 
tion of the implementation). Simulations of low mass planets 
in turbulent discs indicate that the gravitational interaction 
between embedded objects and the turbulent fluctuations 
can have important consequences for the orbital evolution 
jNelson fc PaDaloizoull2004t lNelsonl2005al) . It is for this rea- 
son we computed run N2 to examine the effect this may have 
on the evolution of smaller bodies that also interact with the 
disc gas through the drag force. 




Figure 1. Logarithm of the density in the disc {left panel) and 
vorticity in the local rotating frame (right panel). Both quanti- 
ties are calculated at the end of model Gl. A region of negative 
vorticity is seen in the middle of the disc (at a location r ~ 2.5 
and <f> ~ t/6), well correlated with a density increase on the left 
panel. 



2.2 Model parameters 

The parameters of the models that we computed here are 
described in table We performed four models, labelled 
Gl, G2, Nl and N2. They correspond to increasing particle 
sizes: 5 cm, 25 cm, and 1 m. The first two, for which there is 
good coupling between the gas and the dust, were computed 
with GLOBAL using the two fluid approach. In models Nl 
and N2 the large particles start to decouple from the gas, so 
the N-body approach implemented in NIRVANA was used. 
In the following, we describe the disc model parameters and 
the properties of the MHD turbulence that develops as a 
consequence of the MRI. We then describe the simulation 
results. 



3 RESULTS 
3.1 Disc model 

The turbulent disc structu r e is com puted following the pro- 
cedure described in lNelsonl (|2005a): p varies like r _1 and the 
equation of state is locally isothermal, with c s oc r -0 ' 5 . The 
computational domain ranges in radius from R m in — 1 to 
Rmax = 5 and its vertical extend is AH = 0.28. In order 
to reduce the computational cost of models Gl and G2, the 
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azimuthal domain only covers the interval [0, n/2]. In that 
case, the resolution is (N r , N^,N Z ) = (260, 152, 44). In mod- 
els Nl and N2, which cover the full azimuthal domain, the 
number of grid cells in <j) is four times larger. 

When converting from computational to physical units, 
we assume that r = 2.3 corresponds to 5 AU. The disc mass 
and density are scaled such that the total disc mass would 
be 0.06 Mq within 40 AU (i.e. approximately three times 
the minimum mass nebula). The orbital period at the disc 
inner edge thus corresponds to ~ 6 years. We assume that 
the internal density of solid bodies p B — 3 g cm -3 . 

In all of our models, turbulence is initiated in the disc 
with a non-zero net flux toroidal magnetic field. When the 
MRI starts to saturate, the field is scaled down by a factor 
of two, while t he density profile is reset to its initial value 
jNelsonll2005al) . The simulation is then restarted. This pro- 
cedure gives a typical volume-averaged a value computed 
from the Maxwell and Reynolds stress of ~ 10 -2 . In both 
GLOBAL and NIRVANA, a quasi-steady state is obtained 
at t — 300 (measured throughout this paper in orbits at the 
inner edge of the disc), when the solid component is intro- 
duced. The turbulent state of the disc maintains a rough 
statistical steady state during the run time (~ 180-200 or- 
bits). The disc structure is illustrated by figure Q where 
the left panel shows the density in the midplane of the disc 
at the end of model Gl. A local increase in the density is 
visible in the middle of the disc. To better understand its 
nature, we plot in the right panel of figure^the vorticity to 
which is calculated at each radius in the local rotating frame 
according to: 

w = Vx» 9 . (5) 

Here v g is the difference between the local gas velocity and 
the azimuthally averaged velocity. White areas correspond 
to large positive values of uj and dark zones to large negative 
values. A region with a lower value of uj than the background 
of the disc is seen, which is well correlated with the density 
increase shown on the left panel. This is the s ignature of 
an anticyclonic vortex jjohansen fc KlahJl2005fl . Note that 
the decrease in uj that is created is of smaller amplitude 
than the vorticity fluctuations that result from the MHD 
turbulence. Similar vortices are observed in runs Nl and N2, 
and we observe these structures to survive for the duration 
of the simulations. Their nature and formation mechanism 
are discussed further and compared with previous work in 
section ^] 

3.2 Models Gl and G2 

At time t = 300, the dust is introduced such that p/pd = 100 
everywhere in the disc and the evolution is followed until 
t = 477. The parameter Qt s (where Q is the Keplerian an- 
gular velocity) is independent of radius in these models. Two 
cases are investigated. Model Gl corresponds to Qt s = 0.1, 
and model G2 has f2r s =0.5. Given the model normaliza- 
tion described above, these parameters respectively corre- 
spond to solid bodies whose diameter is equal to 5 and 25 
centimetres. 

During the simulations, the dust drifts radially toward 
the central object and the spatial distribution evolves. The 
logarithm of the gas-to-dust ratio in the (r — <j>) plane is 
shown in figure [5] for model Gl (left panel) and G2 (right 




Figure 2. Logarithm of the dust-to— gas ratio at the end of model 
Gl (left panel) and G2 (right panel). Its maximum value respec- 
tively reaches 0.34 and 1.41. 
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Figure 3. Time history of the maximum value reached by the 
gas to dust ratio for model Gl (solid line) and G2 (dashed line). 



panel) at time t = 176. The striking feature of both of the 
snapshots is the accumulation of dust in a region of the disc 
close to r — 2.7. The peak value of the dust to gas ratio 
at that time is 0.34 and 1.41 for model Gl and G2, which 
indicates an increase of respectively 34 and 141 during the 
course of the simulation. A comparison between figuresQand 
[5] shows that this accumulation occurs at a location corre- 
sponding to a local increase of the gas density (or, equiva- 
lently, of the gas pressure) . We have shown in section l3~T1 that 
this pressure extremum corresponds to an anticyclonic vor- 
tex. The fact that solid bodies are efficiently trapped in such 
a vortex has bee n previously reported in the literature on 
several occasions (Barge & Sommcri al ll995l iGodon fc Livid 
2000; .loha nsen et al.ll2004T) . The novel result of our work is 
that this long-lived vortex is self-consistently generated by 
the MHD turbulence itself. 

Finally, figure [3] plots the time history of the maximum 
value reached on the grid by the gas-to-dust ratio for mod- 
els Gl and G2. Both curves show an increase with time. 
However, there is a clear tendency for larger particles to 
accumulate more quickly in gas pressure maxima, as they 
undergo faster radial migration, and the maximum value for 
the dust-to-gas ratio reached in model G2 is greater than 
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semimajor axis versus time - I metre 
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Figure 4. This figure shows the radial trajectories of 120 rep- 
resentative solid bodies from simulation Nl. Note that there are 
two regions in the disc where particles are trapped, correspond- 
ing to local pressure maxima in the disc. Very similar results were 
obtained in simulation N2. 



unity. In both models, the value of the dust-to-gas ratio 
saturates after about 120 orbits in model Gl and after 50 
orbits in model G2. This is because the outer regions of the 
disc are depleted at that time and there is no more solid 
material to supply the region where dust accumulates. 



One is located at radius r ~ 2.4, and traps 1276 particles. 
Another vortex forms at radius r ~ 2.6, and contains 24 
particles. This vortex is observed to be weaker than the for- 
mer one, explaining the reduced particle concentration. At 
the end of the simulation there are 3 particles orbiting at 
radii r ~ 2.8, but these are not contained within an obvious 
vortex, but appear to be tail-enders whose inward drift has 
been disrupted through turbulence-induced perturbations. 
A third vortex orbits at r ~ 3.0, and contains 341 parti- 
cles at the end of the simulation. Statistically, we find that 
just under 50 percent of the particles migrate substantially 
through the disc, and just over 50 percent become concen- 
trated within the vortices that form. 

A similar picture emerges from run N2, which included 
the gravitational interaction between solid bodies and the 
disc. This simulation resulted in particles becoming concen- 
trated in vortices at radii r ~ 2.5 and r ~ 2.9. The degree of 
concentration was similar to that in run Nl, with 40 percent 
of particles migrating interior to r = 1.5 and 60 percent be- 
ing concentrated further out in the disc. It is clear that gas 
drag dominates the evolution for 1 metre-sized bodies, with 
perturbations to their motion induced by disc gravity having 
little effect. Simulations that include larger bodies (a ^ 10 
metres) indicate that the concentration of solid bodies ob- 
served here does not hold for that size range, and that the 
disc gravity rapidly becomes the dominating influenc e on the 
orbita l evolution of larger planetary building blocks iNelsonl 
12005 til . 



3.3 Models Nl and N2 

Simulations Nl and N2 employed 3000 and 500 particles, 
respectively, to trace the trajectories of 1 metre sized solid 
bodies. Their initial radial and azimuthal locations were dis- 
tributed randomly in a narrow annulus between 2.4 ^ r p 
3.2, and they were initiated with Keplerian circular veloc- 
ities. For the disc model we adopt, the stopping time for 
simulations Nl and N2 are given by Qr s = 1.15 . Be- 

tween radii 1.5 - 3.2, Qt s thus ranges between ~ 1.92 and 
0.85. The radial trajectories of 120 particles from simulation 
Nl are shown in figure |I] The simulation ran for 200 orbits 
at the disc inner edge. The inward drift of particles due to 
gas drag is apparent. For those particles that undergo global 
migration through the disc, we find reasonable overall agree- 
ment in the migration time compared with that obtained in 
an equivalent laminar disc simulation. It is clear, however, 
that there are particles that do not undergo substantial in- 
ward drift over the simulation run time, but instead become 
trapped at specific radii. Figure [I] shows that there are four 
separate radii around which particles become concentrated. 
Most of these particles become trapped because they be- 
come concentrated in anticyclonic vortices that form in the 
turbulent flow, and maintain their identity for the duration 
of the simulation. For a given vortex, the trapped particles 
eventually find themselves being concentrated into a single 
grid cell, such that they then maintain very similar trajec- 
tories for the duration of the simulation. 

Out of 3000 particles considered, 1356 migrate interior 
to r = 1.5 after 200 orbits, without becoming trapped at in- 
termediate radii. There are three distinct vortices that form 
in simulation Nl that lead to particles being concentrated. 



4 DISCUSSION 

We have presented the results of MHD simulations of cylin- 
drical accretion disc models that include the orbital evolu- 
tion of solid bodies ranging in size from 5 centimetres to 
1 metre. A two fluid approach was used for the 5 and 25 
cm sized particles, and an N-body approach used for the 1 
metre sized objects. Both methods were able to follow the 
radial drift of the solid bodies relative to the gas. 

Two results emerge from these simulations: vortices 
were found to develop in the disc and to survive for the 
duration of the simulations, and dust was observed to accu- 
mulate very efficiently in these long-lived structures. 

The first of these findings is significant as previous stud- 
ies of MHD turbulence in cylindrical discs <|Armitagjll998l : 
iHawlevlExnl : ISteinacker fePapaloizoul 12002) " have not re- 
ported long-lived vortices. The reason for this may be that 
these structures are difficult to identify in turbulent discs, 
because they tend to be weak in comparison with the tran- 
sient background features. By-eye inspection of density plots 
often does not reveal vortices, and we found it necessary to 
inspect the vorticity before being convinced that vortices 
were indeed present in our models. 

The issue of the formation and survival of vortices in 
discs is long-standing. Although it has been shown that 
anticyclonic vortices can survive for hundred s of orbits 
jBaree fc SommeriallT99l iGodon fc Lividll99gTl . few mech- 
anisms have been put forward to explain their formation 
in a barotropic fluid. IGodon &i Livid (|2000) suggested that 
small vortices resulting from underlying turbulence could 
merge to form larger, stable anticyclonic vortices. We find 
that turbulence generates small scale variations in vorticity, 
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and it is possible that some of the small scale regions of 
anticyclonic vorticity may merge to form longer lived struc- 
tures. We find that early on in the simulations, when the disc 
is being relaxed prior to the dust component being added, 
larger patches of underlying vorticity also exist. These are 
not apparent against the higher amplitude background fluc- 
tuations, but become so if the magnetic field is switched off 
and the turbulence is allowed to decay. The disc returns to 
a laminar state, but with continued existence of the under- 
lying vortices that are apparently responsible for trapping 
dust when it is added later to the simulations. 

Another possible route to their formation is th e 
Papaloi z ou-Pr ingle instability jPapaloizou fe Pringlell98l . 
lHawlevI (119871) showed that radially-slender tori will break- 
up into gaseous blobs or "planets". There is some evidence 
that the vortices formed in our discs occur near pressure and 
density maxima caused by radial variations in the viscous 
stress. It remains an open question whether these radial den- 
sity variations are susceptible to this instability, causing the 
vortices we observe to form. A detailed analysis of the ori- 
gin of the v ortices goes beyond the sco pe of this letter, and 
we note that lBarranco fc Marcusl d2005T) have suggested that 
vortices are unstable when vertical stratification is included. 
The equatio n of state may also be of som e importance (see, 
for example. iKlahr fc Bodenheimerll2003l) .We will examine 
these issues in a future paper. 

The second result presented in this letter is that 
solid bodies in the size range considered concentrate effi- 
ciently in regions of the disc where pressure maxima form 
due to the anticyclonic vortices described above. This ac- 
cumulation is_a_well_Jcnow_process that has been stud- 
ied before fearge fc Sommerial Il995l : iGodon fc Lrvic1l200ti 
IJohansen et alJ l2004F ~The two fluid calculations demon- 
strated that the effect of trapping was more pronounced for 
the 30 centimetre sized objects than for the 5 centimetre 
ones because of their more rapid inward drift. A simulation 
performed for 1 metre sized bodies including the effects of 
the disc gravity demonstrated that gas drag forces are very 
dominant in determining the evolution for these size ranges. 
All simulations indicated that the trapping was very effi- 
cient, with between 50 - 75 percent of the solids surviving 
inward migration through the full extent of the disc over 
runs times of ~ 200 orbits. Test calculations using laminar 
disc models showed complete loss of solid bodies with sizes 
in the range considered here. 

Note that we did not include the back reaction of the 
solids on the gas in this work, which is valid only in regions 
where the gas density is much larger that the solid density. 
Large enhancements in the density of solids mean that the 
local dynamics may become dominated by the solid bodies 
rather than gas. It will be necessary to include the back re- 
action of the solids on the gas in future studies to investigate 
this effect. 

In the N-body approach , finite resolution effects tend 
to cause particles to clump together on the grid-scale, and it 
is for this reason that we have not quoted density enhance- 
ment factors for these runs. This issue can be addressed by 
adopting a sub grid scale model for the turbulence. 

The effect of trapping presented in this letter may 
help overcome the problem of very rapid migration of 
intermediate sized b odies during planet formation (e.g. 
IWeidenschilliiig||l977h . The enhanced concentration arising 



from this trapping will very likely lead to an enhanced 
growth rate of the solid bodies involved via binary sticking, 
though the sticking presumably occurs between larger bod- 
ies and smaller grains that are trapped in the same region, 
given the difficulty of sticking macroscopic objects together. 
This will allow bodies to grow to sizes for which the inward 
migration is much less rapid. Indeed simulations of larger 
bodies in turbulent discs suggest that objects larger than 10 
metres no longer have their evolution controlled by gas drag 
but rather through gravitational in teraction with turbulent 
density fluctuations (lNelsonll2005ri) . 
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